/*
 * Open Source Physics software is free software as described near the bottom of this code file.
 *
 * For additional information and documentation on Open Source Physics please see:
 * <http://www.opensourcephysics.org/>
 */

package org.opensourcephysics.numerics;
import java.util.ArrayList;

/**
 * Polynomial implements a mathematical polynomial:
 * c[0] + c[1] * x + c[2] * x^2 + ....
 *
 * This class is based on code published in Object-Oriented Implementation of
 * Numerical Methods by Didier H. Besset.
 * The code has been adapted to the OSP framework by Wolfgang Christian.
 */
public class Polynomial implements Function {

  /**
   * Polynomial coefficients.
   */
  protected double[] coefficients;

  /**
   * Constructs a polynomial with the given coefficients.
   * @param coef polynomial coefficients.
   */
  public Polynomial(double[] coef) {
    coefficients = coef;
  }

  /**
   * Gets a clone of the polynomial coefficients c:
   * c[0] + c[1] * x + c[2] * x^2 + ....
   *
   * @return double[]
   */
  public double[] getCoefficients() {
    return(double[]) coefficients.clone();
  }

  /**
   * Constructs a polynomial with the given coefficients.
   * @param coef polynomial coefficients.
   */
  public Polynomial(String[] coef) {
    coefficients = new double[coef.length];
    for(int i = 0, n = coef.length;i<n;i++) {
      try {
        coefficients[i] = Double.parseDouble(coef[i]);
      } catch(NumberFormatException ex) {
        coefficients[i] = 0;
      }
    }
  }

  /**
   * Evaluates a polynomial using the given coefficients.
   *
   * @param x
   * @param coeff the polynomial coefficients.
   */
  public static double evalPolynomial(final double x, final double[] coeff) {
    int n = coeff.length-1;
    double y = coeff[n];
    for(int i = n-1;i>=0;i--) {
      y = coeff[i]+(y*x);
    }
    return y;
  }

  /**
   *
   * @param r double    number added to the polynomial.
   * @return Polynomial
   */
  public Polynomial add(double r) {
    int n = coefficients.length;
    double coef[] = new double[n];
    coef[0] = coefficients[0]+r;
    for(int i = 1;i<n;i++) {
      coef[i] = coefficients[i];
    }
    return new Polynomial(coef);
  }

  /**
   * Adds the given polynomial to this polynomial.
   * @param p Polynomial
   * @return Polynomial
   */
  public Polynomial add(Polynomial p) {
    int n = Math.max(p.degree(), degree())+1;
    double[] coef = new double[n];
    for(int i = 0;i<n;i++) {
      coef[i] = coefficient(i)+p.coefficient(i);
    }
    return new Polynomial(coef);
  }

  /**
   * Gets the coefficient value at the desired position
   * @param n int    the position of the coefficient to be returned
   * @return double the coefficient value
   * @version 1.2
   */
  public double coefficient(int n) {
    return n<coefficients.length ? coefficients[n] : 0;
  }

  /**
   * Deflates the polynomial by removing the root.
   *
   * @param r double     a root of the polynomial (no check made).
   * @return Polynomial the receiver divided by polynomial (x - r).
   */
  public Polynomial deflate(double r) {
    int n = degree();
    double remainder = coefficients[n];
    double[] coef = new double[n];
    for(int k = n-1;k>=0;k--) {
      coef[k] = remainder;
      remainder = remainder*r+coefficients[k];
    }
    return new Polynomial(coef);
  }

  /**
   * Gets the degree of this polynomial function.
   * @return int degree of this polynomial function
   */
  public int degree() {
    return coefficients.length-1;
  }

  /**
   * Gets the derivative of this polynomial.
   * @return Polynomial the derivative.
   */
  public Polynomial derivative() {
    int n = degree();
    if(n==0) {
      double coef[] = {0};
      return new Polynomial(coef);
    }
    double coef[] = new double[n];
    for(int i = 1;i<=n;i++) {
      coef[i-1] = coefficients[i]*i;
    }
    return new Polynomial(coef);
  }

  /**
   * Divides this polynomial by a constant.
   * @param r double
   * @return Polynomial
   */
  public Polynomial divide(double r) {
    return multiply(1/r);
  }

  /**
   * Divides this polynomial by another polynomial.
   *
   * The remainder is dropped.
   *
   * @param p Polynomial
   * @return Polynomial
   */
  public Polynomial divide(Polynomial p) {
    return divideWithRemainder(p)[0];
  }

  /**
   * Divides this polynomial by another polynomial.
   *
   * @param p polynomial
   * @return polynomial array containing the answer and remainder
   */
  public Polynomial[] divideWithRemainder(Polynomial p) {
    Polynomial[] answer = new Polynomial[2];
    int m = degree();
    int n = p.degree();
    if(m<n) {
      double[] q = {0};
      answer[0] = new Polynomial(q);
      answer[1] = p;
      return answer;
    }
    double[] quotient = new double[m-n+1];
    double[] coef = new double[m+1];
    for(int k = 0;k<=m;k++) {
      coef[k] = coefficients[k];
    }
    double norm = 1/p.coefficient(n);
    for(int k = m-n;k>=0;k--) {
      quotient[k] = coef[n+k]*norm;
      for(int j = n+k-1;j>=k;j--) {
        coef[j] -= quotient[k]*p.coefficient(j-k);
      }
    }
    double[] remainder = new double[n];
    for(int k = 0;k<n;k++) {
      remainder[k] = coef[k];
    }
    answer[0] = new Polynomial(quotient);
    answer[1] = new Polynomial(remainder);
    return answer;
  }

  /**
   * Integrates this polynomial.
   * The integral has the value 0 at x = 0.
   *
   * @return Polynomial the integral
   */
  public Polynomial integral() {
    return integral(0);
  }

  /**
   * Integrates this polynomial having the specified value for x = 0.
   * @param value double    value of the integral at x=0
   * @return Polynomial the integral.
   */
  public Polynomial integral(double value) {
    int n = coefficients.length+1;
    double coef[] = new double[n];
    coef[0] = value;
    for(int i = 1;i<n;i++) {
      coef[i] = coefficients[i-1]/i;
    }
    return new Polynomial(coef);
  }

  /**
   * Multipiels this polynomial by a constant.
   * @param r double
   * @return Polynomial
   */
  public Polynomial multiply(double r) {
    int n = coefficients.length;
    double coef[] = new double[n];
    for(int i = 0;i<n;i++) {
      coef[i] = coefficients[i]*r;
    }
    return new Polynomial(coef);
  }

  /**
   * Multiplies this polynomial by another polynomial.
   * @param p Polynomial
   * @return Polynomial
   */
  public Polynomial multiply(Polynomial p) {
    int n = p.degree()+degree();
    double[] coef = new double[n+1];
    for(int i = 0;i<=n;i++) {
      coef[i] = 0;
      for(int k = 0;k<=i;k++) {
        coef[i] += p.coefficient(k)*coefficient(i-k);
      }
    }
    return new Polynomial(coef);
  }

  /**
   * Gets the roots of this polynomial.
   * @return double[]
   */
  public double[] roots() {
    return roots(Util.defaultNumericalPrecision);
  }

  /**
   * Gets the roots of this polynomial.
   * @param desiredPrecision double
   * @return double[]
   */
  public double[] roots(double desiredPrecision) {
    double start = 0;
    if(degree()<1) {
      return new double[0];
    }
    Polynomial dp = derivative();
    int counter = 0;
    while(counter<100&&Math.abs(dp.evaluate(start))<desiredPrecision) {
      start = Math.random();
      counter++;
    }
    Polynomial p = this;
    ArrayList list = new ArrayList(degree());
    while(true) {
      double r = Root.newton(p, dp, start, desiredPrecision);
      if(Double.isNaN(r)) {
        break;
      }
      list.add(new Double(r));
      p = p.deflate(r);
      if(p.degree()==0) {
        break;
      }
      dp = p.derivative();
      start = 0;
      counter = 0;
      while(counter<100&&Math.abs(dp.evaluate(start))<desiredPrecision) {
        start = Math.random();
        counter++;
      }
    }
    double[] roots = new double[list.size()];
    for(int i = 0, n = list.size();i<n;i++) {
      roots[i] = ((Double) list.get(i)).doubleValue();
    }
    return roots;
  }

  /**
   * Substracts a constant from this polynomial.
   *
   * @param r the constant
   * @return Polynomial
   */
  public Polynomial subtract(double r) {
    return add(-r);
  }

  /**
   * Subtracts another polynomial from this polynomial.
   * @return Polynomial
   * @param p Polynomial
   */
  public Polynomial subtract(Polynomial p) {
    int n = Math.max(p.degree(), degree())+1;
    double[] coef = new double[n];
    for(int i = 0;i<n;i++) {
      coef[i] = coefficient(i)-p.coefficient(i);
    }
    return new Polynomial(coef);
  }

  /**
   * Converts this polynomial to a String.
   */
  public String toString() {
    if(coefficients==null||coefficients.length<1) {
      return "Polynomial coefficients are undefined.";
    }
    StringBuffer sb = new StringBuffer();
    boolean firstNonZeroCoefficientPrinted = false;
    for(int n = 0, m = coefficients.length;n<m;n++) {
      if(coefficients[n]!=0) {
        if(firstNonZeroCoefficientPrinted) {
          sb.append(coefficients[n]>0 ? " + " : " ");
        } else {
          firstNonZeroCoefficientPrinted = true;
        }
        if(n==0||coefficients[n]!=1) {
          sb.append(Double.toString(coefficients[n]));
        }
        if(n>0) {
          sb.append(" x^"+n);
        }
      }
    }
    String str = sb.toString();
    if(str.equals("")) {
      return "0";
    } else {
      return str;
    }
  }

  /**
   * Evaluates the polynomial for the specified variable value.
   * @param x double    value at which the polynomial is evaluated
   * @return double polynomial value.
   */
  public double evaluate(double x) {
    int n = coefficients.length;
    double answer = coefficients[--n];
    while(n>0) {
      answer = answer*x+coefficients[--n];
    }
    return answer;
  }

  /**
   * Returns the value and the derivative of this polynomial
   * for the specified variable value in an array of two elements
   *
   * @param x double    value at which the polynomial is evaluated
   * @return double[0]   the value of the polynomial
   * @return double[1]   the derivative of the polynomial
   */
  public double[] valueAndDerivative(double x) {
    int n = coefficients.length;
    double[] answer = new double[2];
    answer[0] = coefficients[--n];
    answer[1] = 0;
    while(n>0) {
      answer[1] = answer[1]*x+answer[0];
      answer[0] = answer[0]*x+coefficients[--n];
    }
    return answer;
  }
}

/*
 * Open Source Physics software is free software; you can redistribute
 * it and/or modify it under the terms of the GNU General Public License (GPL) as
 * published by the Free Software Foundation; either version 2 of the License,
 * or(at your option) any later version.

 * Code that uses any portion of the code in the org.opensourcephysics package
 * or any subpackage (subdirectory) of this package must must also be be released
 * under the GNU GPL license.
 *
 * This software is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
 * or view the license online at http://www.gnu.org/copyleft/gpl.html
 *
 * Copyright (c) 2007  The Open Source Physics project
 *                     http://www.opensourcephysics.org
 */
